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Abstract 



We prove in an abstract setting that standard (continuous) Galerkin finite clement 
approximations are the limit of interior penalty discontinuous Galerkin approxima- 
<**] ' tions as the penalty parameter tends to infinity. We apply this result to equations of 

non-negative characteristic form and the non-linear, time dependent system of incom- 
pressible miscible displacement. Moreover, we investigate varying the penalty parame- 
ter on only a subset of a triangulation and the effects of local super-penalization on the 
stability of the method, resulting in a partly continuous, partly discontinuous method 
in the limit. An iterative automatic procedure is also proposed for the determination 
j of the continuous region of the domain without loss of stability of the method. 

t> 

ir} ! 1 Introduction 

, The discontinuous Galerkin (dG) finite element method has become widely used in recent 

years for a variety of problems as it possesses several desirable qualities, such as: Good 
stability properties due to the natural incorporation of upwinding techniques; flexible mesh 
design as hanging nodes and irregular meshes are admissible; and relatively easy imple- 
^ \ mentation of /ip-adaptive algorithms. These properties however come with the drawback 

of an increased number of degrees of freedom compared to a standard conforming method. 
For instance, when using an axi-parallel quadrilateral mesh in two dimensions with piece- 
wise bilinear elements for which the standard continuous Galerkin (cG) finite element 
method has approximately n degrees of freedom (depending on boundary conditions) the 
dG method on the same mesh has approximately An degrees of freedom. 

For advection-dominated advection-diffusion-reaction equations the standard cG method 
exhibits poor stability properties and non physical oscillations may pollute the approxima- 
tion globally. Discontinuous Galerkin methods have generally better stability properties. 
In the case of interior penalty dG method, for instance, stability in the upwind direction 
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has been shown in the inf-sup sense, e.g., in [TJ [8], generalizing ideas from |18j . where 
purely hyperbolic problems were considered. 

Conceptually, somewhere between the standard cG and interior penalty dG methods 
lies the continuous-discontinuous Galerkin (cdG) finite element method |10| . whereby one 
seeks a Galerkin solution on a finite element space V^, dG with V cG C V cdG C V dG , where 
V cG and V dG are the standard cG and dG finite element spaces. In the context of problems 
with layers or sharp fronts, continuous elements can be used away from the layers/fronts 
and discontinuous elements (accommodating appropriate upwinding) can be used in the 
region where the layers/fronts are present. This idea has been studied previously in the 
context of problems with layers by Dawson and Proft [13] using transmission conditions 
between regions where different spaces are used. Cangiani, Georgoulis and Jensen [10J and 
Devloo, Forti and Gomez [T3] have previously compared the cdG finite element method 
with alternative methods for advection-diffusion equations. 

The control of discontinuities across element interfaces in the dG framework can be ex- 
ercised by introducing and/or tuning the, so-called, jump penalization parameters. Using 
excessive penalization within a dG approximation will be referred to as the super penalty 
method. It is natural to expect that as the penalty parameter is increased the interele- 
ment jumps in the numerical approximation decrease. It has been shown by Larson and 
Niklasson [19J for stationary linear elliptic problems (using the interior penalty method) 
and by Burman, Quarteroni and Stamm [S] for stationary hyperbolic problems (penalising 
the jumps of the approximation for discontinuous elements and the jumps in the gradient 
of the approximation for continuous elements) that the dG approximation converges to 
the cG approximation as the jump penalization parameter tends to infinity. 

In this work, our aim is twofold. Firstly, we present an alternative proof of the con- 
vergence of dG methods to cG methods, using a far more general framework covering the 
cases considered by [£l [19] and also non-linear and time dependent problems. Moreover, 
we show that super-penalization procedures can be localized to designated element faces, 
thereby arriving to partly continuous, partly discontinuous finite element methods. As 
particular examples we consider the limits of the interior penalty dG method for PDEs 
with non-negative characteristic form [17] and the mixed Raviart-Thomas-dG method for 
the miscible displacement system presented in [3]. 

Secondly, we continue the numerical investigations of [10] in the context of blending 
locally continuous and discontinuous methods. In particular, we investigate to what extent 
numerical oscillations appear as local super-penalization is applied. The aim, of course, 
is to find the extent to which degrees of freedom can be removed by using locally con- 
tinuous finite element spaces without affecting the extra stability offered by dG methods. 
To this end, we consider an advection-dominated advection-diffusion problem containing 
boundary layer behaviour, where the continuous and the discontinuous regions of the finite 
element solution are tuned manually. A second example investigates the use of an iterative 
automatic procedure for the determination of the continuous region of the domain by local 
super-penalization without loss of stability of the method. The procedure is applied to 
the problem of incompressible miscible displacement. 

This work is organized as followed. After introducing notation in Section 11.11 an ab- 
stract discussion of the limit of penalty methods is given in Section [21 We then show how 
this framework can be applied to equations of non-negative characteristic form in Section 
[3] and to the non-linear equations of incompressible miscible displacement in Section HJ Fi- 
nally, Section [5] contains a number of numerical experiments and discussion of an iterative 
automatic procedure for determining the continuous regions of the approximation. 
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1.1 Notation 



Let f2 C M. d be a bounded, open polygonal domain. We denote by Th a subdivision of 
into open non-overlapping d-simplices E. The diameter of E G 7^ is denoted by Let 
also £/j := UEeT h dE be the skeleton of the mesh Th, while := £h\dft. Finally, let T 
denote the set of elemental boundary faces, i.e., those which lie in dVL. 

For e G £?, with e = E + (~)E~ for G 7^, we define h e := min(/i B -, h E +). Given 

a generic scalar field v : Q — > R that may be discontinuous across e, we set := v\e±, 
the interior trace on E^ and, similarly, for a generic vector field r : — ► R d . Define the 
average and jump for a generic scalar as 

{M -=\^ + + ^), M ■= v + n + + v-n~, onee4°, 
and for a generic vector field as 

It}} ■= \(t + + O, [t] := r+ • n+ + r" • n" on e G ^, 

where n ± is the outward pointing normal from E 11 * 1 on e. For e G T the definitions become 

:= i/, H := vn, {t} := r, on e G T. 

Given a vector 6 denote the inflow and outflow boundaries of Q by 

d in n = T in :={x G dQ : b ■ n < 0}, 
5 out fi _ pout G ^ : ft . n > o} 

and for an element 

d in E :={x G dE : b ■ n < 0}, 
9°"*^ G d£ : 6-n > 0}. 

We denote the trace of a function v on an edge by v ni (resp. u out ) on the side of the edge 
where b ■ n < (resp. 6 • n > 0). We construct the mesh so that the sign of b ■ n is the 
same for every x G e. 

For the cdG method, we will require the following additional notation. We identify 
a decomposition of our triangulation Th into two disjoint triangulations 7dG and Tg '■= 
T \ TlGj upon which continuous and discontinuous elements will be applied respectively, 
henceforth referred to as the continuous and discontinuous regions of the triangulation. 
Define J := TcG^TdG an d define £dG to be the skeleton of Tig an d £ C G '■= £h\£dG- Note 
that with this definition the faces in J are part of the discontinuous skeleton S^G only. 
Define T c g := Tg l~l T, the set of boundary faces of the continuous region, and similarly 
FdG : = 7dG n r. 

Finally, we will denote by the elementwise divergence operator. 

2 An abstract discussion 

Consider a (possibly non-linear) operator B : W x W — > R where W is a finite dimensional 
vector space with norm || • \\w- Suppose there exists a decomposition of W such that 
V © X = W for V, X C W. In particular this means we can write any w G W uniquely as 
w = v + x for some v G V and x G X. 
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Assume that B is coercive, i.e., there exists A\y > (typically independent of the 
dimension of W), such that 

(2.1) B(w,w) > Aw\\w\\w VwdW. 

Consider another operator S : W x W — > M, whose support is restricted to X x X in 
the sense that 

(2.2) S(v,v) = Vu.ieV 
and 

(2.3) S(y, x) = S(x, v) = V v G V, x G X. 

We require coercivity on X, i.e., there exists Ax > such that for all x G X 

(2.4) S(x,x)>A x \\xf x , 

where is a norm on X. In view of (|2.2p this gives S(w,w) > Ax\\w\\x = Ax||a;||^. 

We construct a further operator 

(2.5) B a :=B + aS 

where < a G K, and call this the super penalised bilinear form. 

Let £ be an element of the dual space W* of W, independent of a. Then choose 
w a G W such that 

(2.6) B a (w a , w) = t(w) V w G W. 
Also choose Vh G V such that 

(2.7) B(t; h ,t;) V v G V. 
Observe that for all <r G M. 

(2.8) B a (vfi, v) = B{v h ,v) = l(v) V « G V 
using Now with (pH]) . flZH) and flZSD we have 

^Wll^vllw/ < B(w a ,w a ) + aS(w a ,w a ) 

= B a (w a ,w a ) 
= £(w a ) 

< \\£\\w*\\w*\\w- 

Using Young's inequality we see 

A 2 1 
(2-9) —\\™*\\w + 2A w A x \\w4 2 x < -\\£\\ 2 W ,. 

a a 

Each of A\y, Ax and ||^||w* are independent of a. We write w a = v a + x a , the unique 
decomposition with v a G V and x CT G A. From (|2.9p we see 

(2.10) lim 1 1 Wo- + x a \x = lim ||x CT ||x = 0. 

Therefore x CT — > as <r — > oo. 



Now assume that B is continuous in the first argument in the following sense: If 
lhrij-j.oo Wi = w G W then 

(2.11) lim B(w»v) =B(w,v) V v G V. 

i— >oo 

Suppose w a v h as a — > oo. Then there exists e > such that there is some sequence 
{ w a(i)}i with a(i) — > oo as i — > oo satisfying 

(2.12) \\w a{t) -v h \\w > e Vi G N. 

Owing to ()2.9p the sequence is a bounded subset of W. Then by the Heine-Borel 

theorem there exists a convergent subsequence, also denoted {w a u\}i, such that 



(2.13) w = lim w a u). 
Considering ()2.10|) we know that w G V. We have that for all v G V 

B(w,v) = B ( lim w c (i),v ) 

= limBK (i))? ;) bydZIU) 

MOO v ' 

= lim B^u;^),*;) by ([23]) 

i — yoo 

= lim £(u) by $ZM) 

i— >oo 

= £(v). 

Hence w satisfies (|2.7p and by ()2.13p we have 

lim \\w a (i) - v h \\ w = 0. 

This contradicts (|2.12p and we conclude that all subsequences {w a (i)}i converge to 
Therefore 

(2.14) lim (w a - v h ) = 0. 

cr— yoo 

We finally remark on the potential loss of stability due to super-penalization. It can 
be seen from (|2.9p that as x a — > when a — > oo the coercivity of £> CT is increasingly 
compromised, which can lead to loss of stability and reduction on the rate of convergence 
in various settings. 

3 Equations of Non-Negative Characteristic Form 

We now examine the diffusion- ad vection-reaction equation (see |17| ) 

^ -V • (A(x)Vu) + b(x) ■ Vu + c(x)u = f(x) in fi, 

u = on d£l 

with b a M. d valued function whose entries are Lipschitz continuous on Q, c G L°°(Q) and 
/ G L 2 (yi) real valued functions. The diffusion coefficient A is a d x d symmetric matrix 
with entries being bounded, piecewise continuous real-valued functions defined on 0, with 

C T AC > VC G M. d , a.e. x G H. 
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Under this condition, (|3.ip is named a partial differential equation with non-negative char- 
acteristic form. 

We define the cdG space to be 

(3.2) V cdG := {v G L 2 (n) : VE G T h , v\ E G P r , v\r cG = 0, v\ TcG G C(T cG )} 

where P r is the space of polynomials of degree at most r supported on E. We define the 
dG space to be 

(3.3) V dG := {v G L 2 (Q) : \/E G T h , v\ E G P r }. 

Finally we define V^ G := V cdG © V ± (corresponding to W := V © X in the notation of 
Section [2]). Note that the standard continuous space is obtained by setting Th = T c g- 

Define B : V dG x V dG — > R, the bilinear form for the interior penalty family of methods 
with i? G {-1, 0, 1} for (|3TTj) . by 

(3.4) u>) := B d (w,w) + B ar (w,w) 
with 

mtw} ■ jwl ds 



(3.5) 



and 



(3.6) 



B d (w,w) := ^2 / AV fe w • V fe -w da; + / , 

Y, J ({{AV^}} • H - ^AV fc t&J • H) ^ 



e&€ h ' 



B ar (w,w) := y^ / (6 • Vhw)w + ewu) da; 
Ber h ^ 



njitwds. 



We define m := C p |Ar 2 }//i e , A := |||\/A| 2 || l°°(e)i with | • 1 2 denoting the matrix- 2- norm, 
and Cp(i9) > fixed for a given $. The linear form is given by 

(3.7) l(w) := f f wdx - 

E£T h JE 

For e G £ c q we have the additional term S : V dG x V dG — > R penalising the jumps where 

(3.8) S(w,w):= [ M H • 1M ds 

{{Ar 2 }}' 



ee£ C G 

and 

M := ( C aT + C7«j 



h e 

with C ar and fixed constants independent of a. Then we define B a (w, w) := B(w, w) + 
aS(w, w). 

Observe that if we take C ar = C d = (or a = 0) we recover the usual interior penalty 
method. If we take C p = C d = C ar = and A = we have the standard (unpenalised) 
bilinear form for the purely hyperbolic equation (assuming of course that we adjust the 
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boundary conditions appropriately). Taking d = and C ar 7^ when A = gives the 
method proposed in [7]. 

All functions in V cdG are continuous on edges in £ c q (recall that by definition edges in 
J are not included in £ c q). Therefore conditions (|2.2j) and (|2.3p are satisfied for this S. 
That is, for any v,v £ V cdG and x G V ± 

(3.9) S{v,v) =S(v,x) =S(x,v) =0. 
We define the following norm for all w G V^ G . 

IMIdG : = E ll^ v ^lli 2 (B) + IIccHI^q) 

(3.10) E6T * , 



where Co := \fc — V^V^ b. We also define for w G V^ G 
(3-11) Hl== E II^MIl! 2 (e)- 

Notice that is a semi-norm on V^ G but a norm on Vj_. To make this distinction clear 
we will write ||a:||s for x G V ± . 

Lemma 3.12. If C p is sufficiently large when # = — 1 then B is coercive on V dG , i.e., for 
all w G V dG 

(3.13) B(w,w)> A cc \\wf dG 
with K\y = 1 when i? = 1 and A^y = 1/2 w/ten •& = —!. 

Proof. See, e.g., [IT] for a proof. □ 

From the definition it is clear that S is coercive with constant one on V ± , i.e., for all 

x G Vj_ 

(3.14) S(x,x) = \\xf s . 

Definition 3.15. Define a dG approximation to (|3.ip as w a G V dG satisfying 
(3.16) £ ff K, w) = %) Vw G V dG . 

Definition 3.17. Define a cdG approximation to (|3.ip as G V dG satisfying 
(3.18) = ^ G V cdG . 

Using (|3.9|) we see that Vh also satisfies B{vh,v) = £(v) for all v G V cdG . 

Theorem 3.19. The dG finite element approximation w a converges to the cdG finite 
element approximation as a — )• oo, i.e., 

lim (w a - v h ) = 0. 

a— >oo 

Proof. Following the argument of Section [5] we use Lemma 13.121 and (|3.14p and note that 
(|2,lip is satisfied as linear operators in finite-dimensional vector spaces are continuous. □ 
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4 Incompressible Miscible Displacement 

We consider the problem of finding the numerical solution to the coupled equations for 
the pressure p = p(t,x), Darcy velocity u = u(t,x) and concentration c = c(t,x) of one 
incompressible fluid in a porous medium being displaced by another. We consider the 
miscible case where both fluids are in the same phase. 

Consider the domain £1? '■= (0, T) x Q. The equations for the miscible displacement 
are given by (e.g., 



(4.1) ip— + u ■ Vc - V • (B(u)Vc) + cq 1 = cq 1 , 

at 

(4.2) V • u = q 1 - q P , 

K 

(4.3) u = —(y.p- p (c)g) 

with the boundary conditions on d£lx '■= (0, T) x d£l given by 

(4.4) u ■ n = 

(4.5) (B(u)Vc) • n = 0, 



and the initial conditions 



(4.6) c(0,-)=c . 

We denote by: <p(x) the porosity of the medium; q 1 > and q p > the pressure at injected 
(source) and production (sink) wells; K(x) the absolute permeability of the medium; //(c) 
the viscosity of the fluid mixture; p(c) the density of the fluid mixture; g the constant vector 
of gravity; D(tt, x) the diffusion-dispersion coefficient; c the injected concentration; and 
Co the initial concentration, which we assume for simplicity to be 0. We define a~ 1 (c) := 
The coupling is nondinear through the coefficients D(u, x), /j,(c) and the advection 
term. We make the common specific choice for the diffusion dispersion tensor, e.g., |12t 

(4.7) D(u, x) = ip (d m I + \u\dtE(u) + \u\d t (I - E(u))) 

where E(u) = uu T /\u\ 2 and I is the identity matrix. We specify that the molecular, 
longitudinal and transverse diffusion coefficients d m , d\ and dt are positive real numbers. 

We solve for the pressure and velocity using a Raviart-Thomas (RT) procedure [151 EO] 
and for the concentration using a cdG method. We refer to the whole scheme as a RT-cdG 
method. For k > we define 

U:={ve (L 2 (n)) 2 : v\ E G (F k (E)) 2 + xF k (E) VE G 7ji, 
v • n continuous on e G 5^}. 

To avoid confusion for the pressure terms we define the space P := V dG where V^ G is 
defined in Then the velocity and pressure are approximated in U x P. To simplify 

the presentation we use the same mesh Th to solve for u, p and c numerically at each time 
step and there is no refinement of the mesh or polynomial degree. However T c g an d 7dG 
are not fixed so the cdG space used to approximate c will vary with time. We define the 
time dependent cdG space by 

Kdc : = i v G l2 ( q ) ■■ yE G T h , v\ E G F r , v\ rj = 0, v\ TJ G C{T{ G )} 
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where 7^ G and r^ G are the T c g region and external boundary of T^ G at time tj. As we 
assert that no change to the shape of the mesh occurs in time we define the time dependent 
dG space as in (|3.3p . Then we may define V^ G := V^ G V[. Note that the de gree k is 
the same for U and P but need not be equal to r, the degree of the polynomials used to 
approximate concentration. 

Let = to < t\ < ••• < ijv = J 1 be a partition of the time interval (0, T). For 
simplicity we assume that each time step is of equal length and define At := tj — t,_i 
and the backward Euler operator D t c^ := (At) _1 (c^ — cj^ 1 ) for j = 1,2, . . . , N. For the 
diffusion part of the concentration equation define the bilinear form 



b*MA) = E o»K)v*i, v fc 4)* + E ™(Ki> k 

(4 8) jBeTh 

- E [(I<iM<)v,4l)e + (K],{{oK)v,<i) c 

for all 4 S V^ dG . The penalty parameter m is defined by [I] 

„ maxjnj B(?x{' + , x)nr, , nj B(n{' , x)ri£, } 
m 2 :£ h ^R, x^ C pcn 1 £h K h ' ± £h V h — h * 

and C pen is chosen such that it is larger than 

sup | h max | itp^ , T^^lr 1 :^6P S ,DG \P*] dxd , E shape regular 



M|| ' \\D^V h v h \\l 



E 



The bilinear form for convection, production and injection is given by the non-standard 
form 

= \ E [(^•Vfc4.4)«-(^V h 4) + ((/ + g p )c(,4) s 
(4-9) 2 ^ h L 

where 4 is defined by 
(4.10) 4'* = 



4 if u{-n + > 0, 
4' + if it{ • n+ < 0. 



This formulation ensures that B cq is semi-definite regardless of the properties of ui. We 
do not need to restrict sums over edges to cells in T c g as i n this region elements of V c &q 
are continuous. Also note that the dG method is a special case of the cdG method where 

r cG - 0. 

Define B(c? h , d J h ;u{) := B d {(P h , d J h ; u{) + B cq {^ h ,d\\ u 3 h ) and 
(4.11) 5(^,4):= [MKl-KJds 



where 



Then for any c^, d° h G V^ dG and x J G V^f 

5(4,4) = 5(c J h ,^)=5( a J,4) = 0. 
We define the following norm for wl G U and G Vj G : 

iikm 2 == E Ha/^(4)v4ii| 2(£) + ^ii«,4ii| a(n) 

(4.12) ^ 

+ E 2\H- n \ 1/2 Knh ( e } + En^KiiiW) 



where go := vr+?' F° r 4 e ^dG define 

(4-13) |4||:= E W^MKMme). 



Notice that (|4.13p is a semi-norm on V^ G but a norm on Vj_. 

Lemma 4.14. If C pcn is chosen large enough then B is coercive for all c? h G V^ G and 
u 3 h G U, i.e., 

(4.15) ^«,<;<)>A W |||4||| 2 . 

Proof. Combine equations (4.3) and (4.6) from [3]. □ 
We have by construction that 5 is coercive with constant one on V7, i.e., for all x J h G 

(4.16) 5(4,4) = ||4|||. 

We discretise the time derivative with the backward Euler operator. Summing over 
each discrete time step gives 

N • N 1 • 1 _! 

E(^>t4>4) = E ^(^4.4) - 1 >4) 

i=i i=i 

AT 

- E ^Il^ 1/2 4llf2(fi) - 2a7 (l^ V2 4 _1 llL2(n) + llv 1/2 4lli»(n)) 

1 / „ y 2 JVi|2 _ II l/ 2 /»°ll 2 



where we have used Young's Inequality. We have assumed that the initial concentration 
is and so Hy^c^llisfn) = 0. 

Definition 4.17. Define the RT-dG approximation (uh,Ph,c a ) G X0 =1 U x HjL-^P x 
H^ =1 V ] dG to (|4.ip - (|4.6p as that generated by the algorithm: For 1 < j < N and c^-" 1 G V^ G 
find (u 3 h) p 3 h ,c' a ) G U x P x V^ G suc/i t/ioi 

(4.18) (V h -« , h ,t4) = (g / -? p ,t4), 

(4-19) (a-v.KVj - (4, • 4) = (K4)5,4) 

/or all (v^jwl) eU x P and 

(4.20) (W,4) +H(4,4;4) + ^(4,4) = («^4) 

/or a// 4 G V^ G . 
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Definition 4.21. Define the RT-cdG approximation (uh,Ph,Ch) 6 n^ 1 C7 x IT^P x 
^-jLiVcdG t° (|4.ip - (|4.6p as that generated by the algorithm: For 1 < j < N and 1 € V^ G 
/mc? (u^p^c^) £U x P x V^ rfG suc/i i/iai 

(4-22) (V fe XXj = (/-(7 P Xj, 

(4-23) (a- x (4)«i.«ft) " (4. V » • 4) = (p(4)<?>4) 

/or a// (Vfrjwjj) £ U x P and 

(4.24) (W„ 4) + 0(4, 4; 4) + ^(4> 4) = («/> rf j 



/or a// 4 G ^ci 



cdG- 



Theorem 4.25. Tfte solution c a E Hj =1 V^ G defined in Definition ^. 17\ converges to Ch € 



TljLiVcdG defined in Definition ^. 21 



j 

as o~ — y oo, i.e., 



(4.26) lim (c a - c h ) = 0. 

G — ^OO 

Proof. Following the argument of Section [2] we use Lemma 14.141 and (|4.16p . In order to 
complete the proof using this argument we must show that for every sequence {c^}j with 
elements in V^ G and lim^oo cj = c 3 G V^ G we have 

(4.27) lim B(4,4;vP(4)) = B(j J^v* {<*)) Vd{ G v[ G 

as in (|2.1ip . where u 3 (-) is the element in U solving (|4.18p - (|4.19|) for a given element of 
V^q. Note that v? : V^ G — > U is a continuous map and so lim^oo u 3 (cj) = u 3 (lim^oo c|) = 
u 3 (c 3 ). This also holds for derivatives as they are taken piecewise. Therefore (|4.2T[) holds 
at each timestep and for the whole discrete solution in time. □ 

5 Numerical Experiments 

We present numerical experiments to illustrate Theorems 13. 191 [4.251 and investigate further 
the performance of the cdG method. 

The results were produced using the CH — \- library deal, ii [21 Ej using both the super 
penalty approach (as a — > oo) and a direct cdG method, i.e., where the test functions are 
in V cdG and therefore by construction there will be no jumps across edges in £ c q. For 
further details of the implementation we refer to [?] . 

5.1 Equations of Non-negative Characteristic Form 

Let 0, = (0, l) 2 . We seek to solve 

-eA« + (l,l) • V« = /. 

Given homogeneous Dirichlet boundary conditions / is chosen such that the solution is 
given by 

/ e (x-l)/ E _ e -l/ E \ / e (y-l)/e _ e -l/e\ 

■=[ X - x-e-i/. [y- x_ e -i/ e • 
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For < e <C 1 this problem exhibits exponential boundary layers along the outflow 
boundaries x = 1 and y = 1 of width 0{e). We consider a uniformly refined mesh of 
squares and set r = 1 (piecewise bilinear polynomials). 

We first look at an example without a layer by setting e = 10. We set T c g = Th, 
i.e., the cG method. Figure [57T1 shows the behaviour of the difference between the dG 
and cG approximations in the L 2 (7~h) norm, ff 1 (7/ l ) semi-norm and the L 2 norm of the 
jumps across edges (represented by [ • J). As a grows the difference in each norm decreases 
linearly. The jumps in the either approximation are already very small, i.e., the dG 
approximation is very close to an element in the cG space. We do not see oscillations 
polluting the continuous approximation. 




a 



Figure 5.1: Example 15.11 with e = 10 and %g — Th- As the penalty parameter is increased the 
difference between the cG and dG approximations decreases linearly in the given norms. 

We now motivate the cdG method by choosing e = 10 -4 and again setting T c g = Th- 
The example now has a sharp layer at the outflow boundaries. We see in Figure I5.2;(a)| 
that increasing a gives a linear response to the error as in Figure 15.11 When we look 
at the error in the dG approximation in Figure I5.2;[b)| we see that the approximation 
becomes worse as the penalty is increased. The layer causes non-physical oscillations to 
pollute the approximation. Although we see convergence of the dG approximation to the 
cG approximation this property is not desirable. 

We now consider the cdG method with e = 5 x 10 -3 and 5 x 10~ 4 , values chosen so 
that the layer is partially resolved in the first case and not resolved in the second. The 
behaviour as a is increased is the same as in the case Th = To, and so we do not plot this. 
We set h = 2 -5 and T c g = (0, 1 — ah) 2 . Varying a£Z determines the number of rows in 
the dG region at the outflow boundary. For e = 5 x 1CP 3 oscillations are apparent in the cG 
approximation but the mesh is sufficiently refined so that they are not large. Decreasing 
a (that is, moving from a fully discontinuous approximation towards a fully continuous 
approximation) results in a small increase in the error of the approximation which can 
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(a) The difference between the cG and dG ap- 
proximations. 



(b) The error in the dG approximation. 



Figure 5.2: Example 15 . 1 1 with e = 10 4 and %g = Th- Now the problem has a layer the error in 
the dG approximation grows as a is increased. Non-physical oscillations pollute the approximation. 

be seen in Figure 15.31 The continuous, discontinuous and continuous-discontinuous ap- 
proximations are very close in the H 1 semi-norm, including when T c g covers the layer. 
When e = 5 x 10 -4 the layer is sufficiently sharp to induce large oscillations in the fully 
continuous approximation. By choosing the T c g region to allow discontinuities in at least 
the final element we may achieve a reduction in degrees of freedom to approximately 30% 
of the discontinuous approximation with very little effect on the error (4096 degrees of 
freedom for the dG approximation, 1024 for the cdG approximation and 1276 for the cdG 
approximation with one row of dG elements at the outflow boundary). 

5.2 Incompressible Miscible Displacement: The "Quarter of Five Spot" 
Problem 

As well as verifying Theorem 14.251 we wish to show that if the region where continuous 
elements are used is chosen appropriately there is little difference in the approximations 
via the RT-cdG or RT-dG method (where the concentration is approximated in the dG 
space). 

We study a standard example [U [TTJ [2T] to illustrate the performance of the cdG 
method for the incompressible miscible displacement problem (|4.ip - (|4.6p . With U = (0, l) 2 
the injection (resp. extraction) well is located at (1,1) (resp. (0,0)). The injection and 
extraction strength are represented over one element by piecewise constant functions such 
that f n q J dx = f n q p dx = 0.018. In g2J we set d x = 1.8 x 10" 4 , d m = 1.8 x 10~ 6 and 
dt = 1.8 x 10 -5 . The porosity is set to 0.1. The concentration dependent viscosity is 
given by //(c) = /^(0)(1 + (.A/f 1 / 4 — l)c) -4 where A4 = 41.0 is the mobility ratio (the ratio 
of the viscosity of the fluids), and fj,(0) = 1. For the initial concentration we set cq = 
corresponding to O uniformly filled with one fluid. Set IK = 0.02881. We consider a uniform 
refinement of Q into squares of side h = 2~ 4 with timestep 4 x 10 -3 and time interval 
(0.0,2.0). With these values a sharp front in the concentration component spreads from 
the injection to extraction point. As can be seen in Figure [5T6^d)| this causes oscillations 
in the continuous approximation. 

First we present the difference between the dG approximation and the cG approxima- 
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Figure 5.3: Example 15 . 1 1 varying T c g for a partially resolved and unresolved layer. In both cases 
the increase in the error as T c g covers the layer (a = 0) is apparent. When the layer is sharper the 
increase is more pronounced as severe oscillations pollute the approximation. 



tion (i.e., with T c g = Th) as a — > oo. In Figure loTH we show \\c a — Ch\\ in both the 1? norm 
against time and the L 2 ((0, T); L 2 (Q)) norm against increasing a. In Figure l5^a)| we see 
a sharp increase in the error over the first few iterations. The initial conditions are in 
the continuous approximation space so the cG and dG approximations are close. As the 
layer spreads through the domain the difference between the cG and dG approximations 
for a given a in the L 2 norm increases slowly. This is because the number of edges in the 
vicinity of the layer increases. Figure [5^b)| shows the same behaviour as the stationary 
examples in Section f5.il 

Picking 7~cG m Section 15.11 was done via knowledge of the true solution and hence 
knowledge of any layers. We do not have this luxury for the problem considered in this 
section. We therefore undertake the following procedure for determining 7^g : 

(1) Determine the initial pressure and velocity given and the injection profile. 

(2) Solve for the first time step using a RT-dG method to find a discontinuous c\. 

(3) For all edges determine || He/,.]] Ilz^a^j. 

(4) Flag every cell where each edge satisfys ||[ch]||x2( e ) < to\. 

(5) If every edge of a cell is nagged set that cell to be part of T c g i n the next iteration. 
Otherwise the element will be in Tag- 

(6) For n iterations use the cdG mesh denned in the previous step. 

(7) For the (re + l) th iteration reset the mesh to be entirely dG, i.e., 7^q +1 = for the 
concentration component, then return to step |(3)[ 
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(a) Evolution of the difference between the cG 
and dG approximations for a = fO 3 to 10 11 in 
L 2 norm. 



(b) Plot of the difference between the cG and 
dG approximations in the L 2 (L 2 ) norm as a is 
increased. 



Figure 5.4: The effect of increasing a for Example 15.21 with T c q = Th- 



The number of iterations between each cdG refinement and the tolerance should con- 
sider the expected motion of the fluid and the time step. We do not consider increasing a 
for the cdG method, but rather study the performance of the method as the tolerance is 
increased by comparing the cdG approximation with a dG approximation where C pen = 10 
and a = 0. With these parameters we set the number of iterations between redefining the 
cdG space to be 5. 

In Figure 15.51 we see that as the tolerance is decreased the difference between the dG 
and cdG approximations in the L 2 norm gets smaller. With a smaller tolerance fewer 
cells are marked as being continuous. The difference introduced by using some continuous 
elements does not seem to propagate in time. 

In Table 15.11 we see that the number of degrees of freedom saved over the simulation 
(500 steps with T = 2.0, At = 4 x 10 -3 ) is considerable. The effect on the approximation 
is however small measured in the L 2 (L 2 ) norm. The number of degrees of freedom for the 
cG method is not 128,000 as would be expected (one degree of freedom per vertex on a 
16 x 16 square mesh for 500 timesteps) due to every fifth iteration being discontinuous. 



tol 


dofs 


\\ C <J - c fe L 2 ((0,T);L 2 (n)) 


cG 


219,470 


3.9970 x 10° 


10~ 3 


323,488 


1.2073 x 10" 2 


10" 4 


355,328 


7.0904 x 10" 4 


10~ 5 


382,384 


1.0455 x 10" 4 


dG 


512,000 


0.0000 x 10° 



Table 5.1: The number of degrees of freedom used for 500 timesteps in Example 15.21 When 
tol = 10~ 5 only 74.6% of the degrees of freedom are used compared to 43% for the continuous 
approximation. 
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Figure 5.5: The behaviour of the cdG approximation for Example I5.2I Using some continuous 
elements does not dramatically increase the error of the cdG approximation compared to the dG 
approximation. 



In Figure loTHl we show the dG, cG and cdG approximations after 380 timesteps. There 
is no visible difference between the plots for dG and cdG at each tolerance (Figures 15.61 



(b) and (c) ). However for the fully continuous approximation the oscillations induced 



by the layer are clearly visible and distort the plot. 
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(c) The cdG approximation with tol = 10 3 . (d) The fully continuous approximation. 

Figure 5.6: A plot of the cG method, cdG method and dG method at time 1.52 (380 time steps). 
The discontinuous region is marked in dark grey for the cdG method. There is no appreciable 
difference between the first three plots. The oscillations are clearly visible in the fully continuous 
plot. 
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